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Resumen 

En este trabajo se present a un metodo de volumenes finitos enriquecido con un 
esquema de multiresolucion completamente adaptativo para obtener adaptatividad es- 
pacial, y un esquema Runge-Kutta-Fehlberg con paso temporal de variacion local para 
obtener adaptatividad temporal, para resolver numericamente las conocidas ecuaciones 
" bidominio" que modelan la actividad electrica del tejido en el miocardio. Se conside- 
ran dos modelos simples para las corrientes de membrana y corrientes ionicas. En 
primer lugar definimos una solucion aproximada y nos referimos a su convergencia 
a la correspondiente solucion debil del problema continuo, obteniendo de este modo 
una demostracion alternativa de que el problema continuo es bien puesto. Luego de 
introducir la tecnica de multiresolucion, se deriva un umbral optimo para descartar la 
informacion no significativa, y tanto la eficiencia como la precision de nuestro metodo 
es vista en terminos de la aceleracion de tiempo de maquina, compresion de memoria 
computacional y errores en diferentes normas. 

1 . Introduccion 

Las mediciones direct as represent an una dificultad obvia en ciencias. Por lo tanto simu- 
laciones numericas son de gran interes, especificamente en modelos cardiacos. Entre tales 
modelos, el modelo bidominio es conocido como uno de los mas precisos y completes para el 
estudio teorico y numerico de la actividad electrica en el tejido cardiaco. Desde el punto de 
vista computacional, el modelo bidominio representa un verdadero desafio, dado que el tejido 
cardiaco tiene tamanos del orden de centimetres y por ejemplo, los frentes de exitacion de las 
ondas son del orden de los 10~^ centimetres. Esta caracteristica local, no solo espacial, sino 
tambien temporal, junto con la aparicion de frentes perfilados en el campo de los potenciales 
elect ricos, hace que las simulaciones en mallas uniformes sean practicamente imposibles de 




Uevar a cabo. Aqui es donde los metodos adaptativos juegan un rol vital en las simulacio- 
nes cardiacas. En este articulo desarroUamos un esquema de multiresolucion completamente 
adaptativo provisto de adaptatividad temporal a traves de una estrategia de paso temporal 
local y deducimos un umbral optimo para descartar informacion no significativa. Basados en 
experiencia previa sobre sistemas de reaccion-difusion y ecuaciones parabolicas degeneradas 
[2, 4, 5], sugerimos que la multiresolucion puede ser una herramienta efectiva para resolver 
las ecuaciones del modelo bidominio. En el contexto de metodos de multiresolucion comple- 
tamente adaptativos, mencionamos que existen trabajos desarroUados por varios grupos de 
investigacion (ver [2, 4, 6, 10, 11]), con aplicaciones enfocadas a otras areas. 

1.1. El modelo bidominio 

Supongamos J"^ C abierto y acotado con frontera suave dQ. En este modelo, el musculo 
cardiaco es interpretado como la union de dos medios continuos interpenetrados y superim- 
puestos: el medio intracelular y el extracelular. Estos ocupan la misma area y se encuentran 
separados por la membrana celular cardiaca. Ui — Ui{t^x) y Uq — UQ{t^x) representan los 
potenciales electricos intracelular y extracelular en G := x (0,T), y la diferencia 
entre estos potenciales v = v{t,x) = conocido como el potencial transmembrana. 

La conductividad del tejido est a represent ada por tensores escalados Mi(x) y Me(x) dados 
por 

Mj(x) = a^I + {al - a^)^i{x)^f{x), 

donde a- = crj (x) G C^(R^) y = crj (x) G (7^(]R^), para j = i, e son las conductividades intra- 
y extracelulares a lo largo y a traves respectivamente de la direccion de la fibra muscular 
correspondiente (paralela a a/(x)). 

Comunmente se utilizan los radios de anisotropia ^ y ^. En general las conductivida- 
des en la direccion longitudinal / son de mayor magnitud que aquellas a traves de la fibra 
(direccion t); y tal caso se denomina anisotropia fuerte en conductividad electrica. 

El siguiente sistema fuertemente acoplado de reaccion-difusion forma el modelo bidominio 
anisotropico (ver [12]): 

(3CmdtV + V • (Me(x)V2ie) + /5/ion(^, w) = ^pp CU f^r, 

-V • ((Mi(x) + Me(x))Vi/e) - V • (Mi(x)Vi;) = en (1) 

dfW — H{y^ w) —{) en Qt- 

Aqui, Cm > representa la capacitancia de superficie de la membrana, es la razon area- 
volumen, y w{t^x) es la variable de recuperacion, que toma en cuenta las variables de con- 
centracion del modelo. Las corrientes de los estimulos aplicados a los medios intra- y extra- 
celulares estan representadas por la funcion /app = Iapp{t^ x) que satisface Iapp{t^ x) dx — 
para casi todo t G (0, T). La funcion H en la ecuacion diferencial ordinaria de (1) y la funcion 
/ion corresponden a uno de los modelos mas simples para las corrientes de la membrana y 
corriente ionica (entre una amplia variedad de tales modelos): el modelo de membrana de 
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Mitchell-Schaeffer [9] 



H{v^ w) 



Woo{v/Vp) - W 




VpV2 
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V 



'^(1 — v/vp)w 



donde 



r/3 para 5 < 7/5, 
r]4 en otro caso 



1 para 5 < 775, 
en otro caso. 



Rm representa la resistividad superficial de la membrana y '^p, /yi, 772, ^3, ^4, ^5 son constantes 
dadas. El sistema (1) es provisto con condiciones de borde de no-flujo, represent ando un 
tejido cardiaco aislado 



y condiciones iniciales apropiadas en Q para el potencial transmembrana y variable de re- 
cuperacion v{0^x) = Vo{x)^ w{0^x) = Wo{x). Para asegurar la dependencia continua de los 
datos en la component e v de la solucion, requerimos que el dato inicial Vq sea compatible 
con (2). Por lo tanto la condicion de compatibilidad 



debe ser satisfecha. 

La teoria estandar para ecuaciones parabolicas-elipticas no puede ser aplicada de forma 
natural en el analisis de las ecuaciones del modelo bidominio, debido a la diferencia existente 
entre los grados de anisotropia entre los medios intra- y extra celulares. Debido a esta carac- 
teristica, el sistema (1) es de naturaleza parabolica degenerada. En [3] los autores prueban 
existencia y unicidad de solucion para las ecuaciones de bidominio, utilizando el metodo de 
Faedo-Galerkin y teoria de compacidad. 

2. Un metodo base de volumenes finitos 

Para resolver numericamente (1) introducimos un metodo estandar de volumenes finitos. 
Una malla admisible para Q sera formada por una familia T de volumenes de control (poligo- 
nos abiertos y convexos) de diametro maximo h. Para todo K ^ T ^ xk denota el centro de 

N{K) el conjunto de vecinos de £int{K) el conjunto de bordes de K en el interior de 
^ y ^ext(^) el conjunto de bordes de K sobre la frontera dQ. Para todo L G N{K) d{K^ L) 
denota la distancia entre xk J Xl^ (Jk,l es la interfaz entre K y y r]K,L {VK,a respect iva- 
mente) es el vector unitario normal a aK,L orientado desde K hacia L. Para todo K ^ \K\ 
es la medida de K. La admisibilidad de T implica que Q = UxerK^ KDL — ^ si K^L^T 
y K ^ y ademas existe una sucesion finita {xK)KeT^ tal que xkXl es ortogonal a (Jk^l- 
Ahora, considerar X G T y L G N{K) con vertices comunes (a^,K,L)i<^</ con / G N\{0}, 
y denotemos por Tk^l al poligono abierto y convexo de vertices (xk^Xl) y {ai^K,L)i<i<i- 
Sea V una discretizacion admisible de Qt, que consiste en una malla admisible para Jl, un 



(Mj(x)Vuj) • n = sobre Sr := d^} x (0,r), j = i,e, 



(2) 




(3) 
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paso temporal At > 0, y N > elegido como el menor entero tal que NAt > T. Con esto, 
escribimos = nAt para n G {0, . . . , A^}. Sobre cada elemento K E T, se definen tensores 
de conductividad (definidos positives) mediante 



J Mj(x)rfx, j = i,e. 



Sea Fj^K,L una aproximacion de J^^^ Mj(a;)VMj • r]K,Ldj para j = i, e, y para K e Qr, 
L e N{K), sea 



■■= \Mj,K Vk,l\ , j = i,e. 



Los flujos difusivos Mj{x)Wuj • r]K,L sobre (7k,l son aproximados por 

XI' 

(Mj(a:)V%) • riK,Ld-i « \(JK,LWi,K,L-f. 



donde y^j es el centro de gk^l Y es una aproximacion de u^{ya)^ ] = h^- La conservatividad 
del metodo nos permite determinar las incognitas adicionales u^^^j^y ademas calcular los flujos 
numericos sobre los hordes: 



Fj,K,L = ^Ik.l ^^j^\^ - ^j,^) si L G N{K), 

donde 



dlK,L = ,,r. - ' Mr\ diK, L). 



d{K, aK,L)M;^K,L + d{L, aK,L)M;^L,K 
Finalmente, aproximaremos el sistema (1) mediante la siguiente formulacion de volumenes 
finitos: Determinar {ulj^)KeT para j = i, e y n G {0, . . . , A^}, {v]^)KeT = (K^k - K,K)KeT 
para n e {0, ...,A^}, y {w^)KeT para n G {0, ...,A^}, tales que para todo K e T y 
n e {0,...,7V- 1} 

^K = jj^^ M^) dx, ^K = jj^J^ wo{x) dx, (4) 

PCJK\ ^ ^ + E <i.,Lj^KL-<,K)+/?l^|/ron,i. = l^|/rpp,K, 



LeN{K) 

\^k,l\ 



LeN(K) ^ ' ^ 



At 



La condicion de horde (2) es tomada en cuenta imponiendo condiciones de no-flujo sohre 
los hordes externos: 

dlK,a-^0^^KL - <k) = for a G Se^K), j = i, e, (5) 
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y (3) es discretizada mediante "^^er l^l'^e ^ P^^^ ^^^^ nG{l,...,A^— 1}. 

La existencia, unicidad de solucion aproximada y convergencia del esquema numerico 
hacia la solucion debil correspondiente, son analizados en el trabajo [1]. Aiin mas, como en 
[2], es posible deducir que en el caso de mallas uniformes, el esquema (4)-(5) es estable bajo 
la condicion CFL 

A . . A ^ 1 A — A- 

2max(^|/ion,K| + 2|4pp,k| j + 4/i max(^|M^,K| + \Me,K\) 

3. Multiresolucion y wavelets 

Considerar como dominio computacional, un simple rectangulo que luego de un cambio 
de variables corresponde a ^2 = [0, 1]^. En primer lugar introducimos una jerarquia de mallas 
anidadas Aq C • • • C A^,, usando una particion diadica uniforme de ft. Cada malla A/ := 
{V{ij),i}{ij)^ con (i^j) a ser definido, esta formada por volumenes de control en cada nivel 
de resolucion V^ij)^i := 2-^[i,i + 1] x [jj + 1], i, j G // = {0, . . . , 2^ - 1}, / = 0, . . . , L. / = 
corresponde al nivel mas grueso y / = L al mas fino. Introducimos tambien los conjuntos de 
refinamiento M(ij)^i = {2(i, j) + e}, e G £^ := {0, 1}^ con :^M(ij)^i = 4. Para cada nivel 
/ = 0, . . . , L, se define la funcion de escala 

= |..^ | X%,),,(x) = 2^^X[o,i]2(2^^i -i,2^X2 - j), 

y por lo tanto el promedio de u{-^t) G L^{Q) sobre el volumen de control V(^,j),; puede 
expresarse en terminos del producto interior 'U(^j)^/ := (t^, (^(i,j),/)Li(i7)- Con esto en mente, 
es posible definir una relacion de dos escalas para las funciones de escala y medias en celda 
respectivamente 

Tal relacion define un operador de proyeccion que transforma elementos desde niveles finos 
a niveles gruesos. Para x G V2(z,j)+a,/+i5 a G definimos la funcion wavelet en funcion de la 
funcion de escala sobre un nivel mas fino, como 



W(iJ),e,l = (-1) (/:?2(i,i)+a,/+l = 2^ \ (-1) ¥^r,/+l. 

Por otro lado, para todo e G := E\{{0, 0)}, es posible obtener una relacion de dos escalas 
inversa (ver [10]) 

Ahora, para e G introducimos los detalles^ que juegan un papel crucial en la deteccion 
de zonas donde la solucion posee altos gradientes 
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Dada la relacion de dos escalas inversa, es posible escribir: 



= 9[iJ),rUr,h (7) 

donde S[.^j^ •={^([V2]+ri,[i/2]+r2),/}ri,r2G{-.,...,o,...,.} ^enota el estencil de interpolacion, ^[.^.^^ 
son coeficientes de la interpolacion, y el gorro sobre u en el lado izquierdo de (7) indica que se 
trata de un valor predicho. La relacion (7) define un operador de prediccion, que transforma 
elementos de niveles gruesos en niveles finos. En contraste con el operador de proyeccion, el 
operador de prediccion no es unico. Sin embargo, se imponen dos restricciones basicas para 
su definicion: Debe ser consistente con la proyeccion, en el sentido de que es el inverso a 
derecha del operador de proyeccion; y debe ser local, en el sentido de que el valor predicho 
dependera solo de j)- En nuestro caso particular, la prediccion es polinomial: 

donde ei, 62 G {0, 1} y 

s s 

Qx • ^ ^^ ^n{^i-\-n,j,l "^1—71^,1^ Qy • ^ ^^ ^pi^i,j-\-p,l '^i,j—p,l\ 
n=l p=l 
s s 

Qxy • ^ ^ ^ ^ ^p{^i-\-n^j-\-p^l '^i-\-n,j—p,l '^i—n,j-\-p,l ~l~ "^i—n^j—p^l)- 
n=l p=l 

Los coeficientes correspondientes son 71 = Y 72 — jf^- Cuanto mas regular es la 

funcion u sobre V(ij)^/, el coeficiente de detalle correspondiente es mas pequeno en modulo. 
En vista de esta propiedad de cancelacion, es natural pensar en alguna estrategia para 
eliminar informacion no significativa {estrategia de corte). La idea basica es eliminar todos 
los elementos de la malla que correspondan a detalles que se encuentran bajo una tolerancia 
(dependiente del nivel de resolucion) dada por 

ei = 2^('-^)s^, (8) 

donde es una tolerancia de referenda a ser determinada en la seccion 3.2. 



3.1. Estructura de datos en arbol 

Organizaremos las medias en celda y los detalles correspondientes utilizando una estruc- 
tura de arbol graduado dinamico. Este tipo de almacenamiento garantiza la estabilidad de 
las operaciones multiescala (ver [6]). Llamaremos raiz a la base, y nodo a cada elemento del 
arbol. Un nodo padre posee cuatro hijos, y un nodo sin hijos es Uamado hoja. Cada nodo 
posee s' — 2 vecinos en cada direccion espacial, Uamados primos cercanos^ necesarios para 
determinar los fiujos en cada hoja; si tales primos cercanos no existieran, deben ser creados 
artificialmente como hojas virtuales. Las hojas del arbol son los elementos que conforman la 
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malla adapt ativa. Denotamos por C{A) a la restriccion del conjunto de nodos A al conjunto 
de hojas. 

En cada paso temporal, la reconstruccion de multiresolucion es aplicada a la parte espacial 
de la solucion u = (v^Uq^w). Luego de cada paso temporal, es necesario actualizar la malla 
adaptativa, y esto se realiza mediante la aplicacion de la estrategia de corte. Una vez que se 
aplica tal estrategia, es necesario agregar una zona de seguridad a la estructura de arbol con 
el fin de asegurar que la malla a ser utilizada en el paso temporal siguiente, sera adecuada 
para represent ar la solucion correspondiente. Est a zona de seguridad sera implement ada 
agregando un nivel mas fino en todas las posiciones posible (sin destruir la estructura de 
arbol). Para cuantificar la mejora obtenida en compresion de datos y tiempo de maquina, 
usaremos la tasa de compresion de datos y la aceleracion de tiempo de maquina (ver [5]) 

_ Af _ CPUtimepv 

2-(^+W + #/:(A)' CPUtimcMR* 

Aqui JV es el numero de element os en la malla mas fina en el nivel L, y #£(A) es la 
cardinalidad del conjunto de hojas. 



3.2. Analisis de error para el metodo de multiresolucion 

Usando las propiedades basicas del esquema de volumenes finitos de referenda, derivamos 
la eleccion optima para la tolerancia de referenda (8). En primer lugar, el error global entre la 
solucion de referenda y la solucion mediante multiresolucion es descompuesto en dos errores 

^ex ~ "^MR 11 — 11 ^ex ~ "^FV 1 1 + 1 1 ^FV ~ "^MR 1 1 ' 

El primer error del lado derecho es denominado error de discretizacion y el segundo error 
es denominado error de perturbacion. Utilizando estimaciones estandar para ambos errores 
y la condicion CFL (6), obtenemos (ver detalles en [4, 6]) que si la tolerancia de referenda 
es dada por 

2(2-a)L-2 

max(^|/ion,K| + 2|/app,K|) + D m^x(^\Mi^K\ + |^e,x|) 

entonces el error de discretizacion y el error de perturbacion poseen el mismo orden de 
magnitud. 



3.3. Aceleracion de la evolucion temporal 

La idea basica del metodo a ser presentado, es utilizar una condicion CFL local, impo- 
niendo el mismo numero CFL para todas las escalas. La estrategia consiste en evolucionar 
todas la hojas situadas en el nivel / usando el paso temporal local 

Ati^2^-^At, / = L-1,...,0, 
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Figura 1: Potencial transmemhrana v, potencial extracelular Uq, y hojas de la correspondiente 
estructura de drbol en los instantes t = 0,01yt=l,l. 



donde At = AIl corresponde al paso temporal sobre el nivel mas fino L. Tal estrategia 
permite incremental el paso temporal para la mayor parte de la malla adaptativa, sin ne- 
cesidad de violar la condicion CFL. Ahora, es necesario sincronizar el paso temporal para 
las porciones de la solucion que se encuentran en distintos niveles de resolucion. Pero esta 
sincronizacion es alcanzada de manera automatica al cabo de 2^ pasos usando At/. 

La proyeccion y prediccion de multiresolucion son efectuadas solo en los niveles ocupados 
por hojas del arbol correspondiente, y solo cada dos pasos temporales. Para el resto de los 
pasos intermedios, se utiliza la misma estructura de arbol. Del mismo modo, los flujos son 
calculados solo en niveles que contienen hojas, y estos son calculados del siguiente modo: 
Si un borde dado es compartido por dos hojas en el mismo nivel /, entonces el calculo del 
flujo se hace de manera estandar, utilizando los primos cercanos u hojas virtuales si fuere 
necesario. Si el borde es compartido por una hoja en el nivel / y una hoja en un nivel mas 
fino / + 1 {borde interfaz)^ calculamos los flujos en el nivel / + 1 sobre el mismo borde, y el 
flujo correspondiente en el nivel / sera igual a la suma de los flujos (en la direccion opuesta) 
sobre los hijos correspondientes en el nivel / + 1. Con el fln de tener siempre a disposicion los 
flujos calculados, la estrategia de paso temporal local debe realizarse recursivamente desde 
el nivel mas flno hasta el nivel mas grueso. 
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Cuadro 1: Tiempo de simulacion, aceleracion de tiempo de mdquina V , tasa de compresion 77 
y errores normalizados. 

4. Ejemplo numerico 

En las simulaciones se utiliza un dominio computacional simple Q = [0, 1 cm]^ y los 
siguientes parametros (siguiendo [12]): capacitancia de la membrana Cm = 1,0 mF/cm^, con- 
ductividades a? = 6i7~^cm"^, al = 0,6 J7~^cm~^, = 24f^~^cm"^ y al = 12 fi~^cm~-^, 
razon superficie volumen P = 4036,5 cm~^, resistividad de superficie Rm = 2 x lO^ricm^, 
Vp — 100 mV, 771 = 0,005, 772 = 0,1, 773 = 1,5, 774 = 7,5, 775 = 0,1. Las fibras forman un angulo 
de — 7r/4 con el eje x y como dato inicial, aplicamos un estimulo en el medio extracelular en 
el centro del dominio (ver Figura 1). Se elige la siguiente configuracion para el metodo de 
multiresolucion: Wavelets con r = 3 moment os nulos, nivel maximal de resolucion L = 9 y 
por lo tanto una malla fina de Af — 65536 elementos, una tolerancia de referencia dada por 
Sr — 5,0 X 10~^. Mostramos en las Figuras 1,2 una secuencia de instantaneas de la solucion 
despues de haber aplicado un estimulo en el centro del dominio. Se muestra tanto la solucion, 
como la correspondiente malla adaptativa generada por la multiresolucion. Los errores han 
sido calculados utilizando como referencia, una solucion aproximada de voliimenes finitos 
sobre una malla con Af = 1024^ = 1048576 voliimenes de control. En la tabla 1 puede notar- 
se que la solucion numerica obtenida aplicando multiresolucion es suficientemente precisa 
(errores del orden de 10~^) y las tasas de compresion son considerablemente altas. 

Para la integracion temporal usando LTS, elegimos CFLq = 0,5 para el nivel mas grueso 
de resolucion, y CFL/ = 2^CFLo para los niveles mas finos. Utilizando LTS, se obtiene un 
aumento sustancial en tasa de aceleracion, y sin embargo los errores se mantienen con el 
mismo orden de precision. 
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Figura 2: Potencial transmembrana v, potencial extracelular Uq, y hojas de la correspondiente 
estructura de drbol en los instantes t = 2^2 y t = 3,3. 
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